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We show that the visibility in interference experiments with Bose-Einstein condensates is directly 
related to the condensate fraction. The probability distribution of the contrast over many runs of 
an interference experiment thus gives the full counting statistics of the condensed atom number. 
For two-dimensional Bose gases, we discuss the universal behavior of the probability distribution in 
the superfluid regime and provide analytical expressions for the distributions for both homogeneous 
and harmonically trapped samples. They are non-Gaussian and unimodal with a variance that is 
directly related to the superfluid density. In general, the visibility is a self- averaging observable only 
in the presence of long range phase coherence. Close to the transition temperature, the visibility 
distribution reflects the universal order parameter distribution in the vicinity of the critical point. 
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I. INTRODUCTION 

Interference experiments constitute an invaluable tool 
for the characterization of the coherence properties of ul- 
tracold gases [TH3]. These properties are particularly in- 
triguing in the case of ultracold one- or two-dimensional 
Bose gases @rE]. Due to strong phase fluctuations, both 
the Id gas at zero temperature and the 2d gas at finite 
temperature exhibit only quasi- long-range order, i.e., 
the one-body density matrix which measures long range 
phase coherence decays as a power law instead of con- 
verging to a finite constant for long distances. Due to 
the quantum nature of the interfering matter fields, the 
measured visibility is a random variable that differs from 
one experimental run to the other. For the case of one- 
dimensional Bose gases at zero temperature, the proba- 
bility distribution of the interference contrast has been 
calculated analytically for arbitrary strong interactions, 
using a mapping to the exactly solvable boundary sine- 
Gordon theory [71. In the weak coupling limit, it turns 
out to be a Gumbel distribution [5J |S] . Numerical data 
calculated within the same theoretical framework, but at 
finite temperature, are in good agreement with experi- 
ment [TU]. 

Here we reconsider the issue of the statistics of the in- 
terference contrast for Bose-Einstein condensates in ar- 
bitrary dimension, discussing in particular the 2d case 
and the connection between long range order and self- 
averaging. We show that the probability distribution 
of the interference contrast is identical to that of the 
condensate fraction in the limit of a large integration 
volume in the absorption images. The resulting distri- 
butions therefore provide the precise counting statistics 
for the number of condensed atoms [TT] . The statistics 
of the condensate number has a universal form in two 
limiting cases: at the critical point as a consequence 
of finite size scaling [T^l [T3] and at temperatures far 
below the critical temperature. We discuss both cases 
and give analytical expressions for the distribution for 
the latter. Specifically, for 2d Bose gases, the distribu- 



tion at temperatures T far below the critical tempera- 
ture of the Berezinskii-Kosterlitz-Thouless (BKT) tran- 
sition [TU HI] , is controlled by the dimensionless param- 
eter T)(T) = l/n s Ay, where n s is the superfluid density, 
and Xt = \J 2nh 2 / mk-QT the thermal wavelength. The 
fluctuations around the average visibility are determined 
by the superfluid density, a quantity that is rather diffi- 
cult to measure by other means [16] . For a homogeneous 
square sample, we show that the probability distribution 
of the interference contrast is close to a convolution of 
two Gumbel distributions, similar to but different from 
the Gumbel distribution that is obtained for the weakly 
interacting limit of a Id Bose gas at zero temperature 
[SI E] - Non-Gaussian distributions are also found in har- 
monically trapped and strongly anisotropic 2d gases, in 
qualitative agreement with preliminary data taken at the 
Ecole normale superieure (ENS) in Paris j!7j . The prin- 
cipal focus of this article is on the physics of 2d Bose 
gases, but we will highlight the differences and similar- 
ities to the 3d case as well as the Id case at vanishing 
temperature. 

This article is organized as follows. In section [TTJ we 
introduce the physical system and discuss the connec- 
tion between the measured distribution of the interfer- 
ence contrast and that of the condensed fraction. In 



section III we use a functional integral description for 
explicit analytical or numerical calculations of the visi- 
bility distribution in the regime where phase fluctuations 
are dominant. In particular, we discuss the general form 
of the probability distribution in terms of its cumulants 
and the issue of self-averaging. Section |IV| is devoted 
to the explicit analytical calculation of the probability 
distribution in the 2d case far below the BKT transi- 
tion temperature for both a homogeneous system and a 
harmonically trapped sample. In section |Vj we discuss 
the scaling behavior of the probability distribution of the 
interference contrast at the critical point, where the av- 
erage visibility vanishes. The one dimensional case and 
the anomalous fluctuations of the condensate fraction in 
three dimensions are discussed in the appendices. 
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FIG. 1: (Color online) Illustration of the experimental situ- 
ation in the 2d case. After a time t of rapid expansion, two 
identical Bose gases overlap and interfere. The interference 
pattern is characterized by its iso-phase surfaces. Three sur- 
faces having all the same phase are shown for illustration, 
the pattern itself continues all over the vertical extension of 
the sample, as suggested by the continuation dots above and 
below the shown surfaces. The sample is imaged using an 
absorption beam in the y direction. 



II. INTERFERENCE STATISTICS 

Typical interference experiments with 2d Bose gases 
{e.g., 0HH]) start out by preparing a pair of such gases 
confined to the lowest transverse mode in the z direction 
and separated by a distance d in the z direction. The 
atoms are then released from the trap and imaged after 
an adjustable free expansion time using absorption imag- 
ing (time of flight measurement). When the trapping 
potentials are cut off, the gases rapidly expand in the 
z direction while the density distribution as a function 
of x and y can approximately be regarded as constant. 
Within this initial expansion period, there is no trans- 
formation of phase fluctuations into density fluctuations, 
which only sets in at a later stage |19j . Due to the rapid 
expansion along z, the gases completely overlap after a 
time t of the order of 20 ms for typical traps. The dif- 
ference between their phases then results in a spatially 
modulated interference pattern (see Fig. [TJ. 

The operator for the atomic density after time of flight 
at an observation point x can be written as [20] 



h{x) = n {z) hi{r) + n 2 {r) + A{r)e lQz + h.c. 



(1) 



Here, tiq{z) is an envelope function {i.e., a normalized 
function of z which has a negligible Fourier component 
at the wave number Q of the expansion), hi 2 are f ne 
in situ {i.e., before time of flight) density operators of 
the individual gases, Q = md/ht the wave vector associ- 
ated with the ballistic expansion and A{r) = ^{(r)^^) 
the operator that determines the local interference am- 
plitude. Here and in the following, we use r to denote a 
point in the trap before time of flight, and x for obser- 
vation points after time of flight [BJ. 



In a given experimental run, the measured density dis- 
tribution n{x) corresponds to an eigenstate of the her- 
mitian operator h{x) and AM can be replaced by a 
complex number A = V{r)e l ^ r \ The resulting den- 
sity can then be written in the standard form no{x){l + 
V{r) cos{Q[z — z {r)]}) of an interference pattern with 
the local visibility V{r) — \A\ and a spatially varying 
shift zo{r) = —<fi{r)/Q. When the interference pattern is 
integrated over a finite volume, spatial variations of both 
V{r) (caused by density fluctuations) and zo{r) (caused 
by phase fluctuations) lead to a reduction of the inte- 
grated visibility V. The height function zq defines a sur- 
face in real space, the iso-phase surface. In one or two 
dimensions, there is one unique iso-phase surface which 
is repeated over the entire extension of the sample (cf. 
Fig. [T]) . In three dimensions, zq depends on all three spa- 
tial coordinates and the shape of the iso-phase surfaces is 
different for different phase values. Provided that density 
fluctuations are negligible, which is always the case in the 
strongly degenerate regime [21], the integrated interfer- 
ence contrast is completely determined by the shape of 
the iso-phase surfaces. 

The interference amplitude in a given run of the exper- 
iment can be extracted from the measured density n{x) 
by taking the Fourier transform along z and evaluating 
it at the wave vector Q, where the magnitude of the local 
interference amplitude 



Aq{x,v) = / dzn(x)e 



iQz 



(2) 



takes its maximum. This yields a complex number which 
contains the random relative phase between the two 
clouds. Its average over many runs will therefore van- 
ish. Here we are interested in the modulus square of 
Aq, which determines the observed visibility of the in- 
terference fringes. Experimentally, the local amplitude 
\Aq{x, y)\ 2 is not a directly accessible quantity since the 
absorption imaging automatically integrates over the y 
direction. In practice, the averaged interference contrast 
is obtained by extending the domain of integration over 
a volume Q that typically covers the entire sample. It is 
then convenient to define an operator 

&= [ d 3 x [ dV h{x)h{x')e- lQ{z - z,) , (3) 
Jn Jn 

whose eigenvalues represent the magnitude of the inte- 
grated contrast in an individual run |20) . In terms of 
the basic in-trap field operators, this operator can be ex- 
pressed by 



d d rd a 



r' ^l{r)Mr')MrW 2 {r') 



(4) 



Now, in a homogeneous condensate, the one-body den- 
sity matrix tpl{r)ipi{r') approaches a constant on length 
scales \r — r'\ larger than the healing length. Moreover, 
one has ^{t^^t') ~ Y2i r ')'4 ) 2[ r ) up to corrections that 
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vanish like l/f2. As a result, the operator 



d~ fT 2 f / d d rdV $(r)<0i(r') 



x d d rd d r'^ 2 (r)Mr')j = N^N^ , (5) 

whose eigenvalues determine the measured interference 
contrast in a given run, is equal to the product of the 
number of condensed atoms 



N, 



(0 



Q- 



dVdV ipl(r)il)i(r') 



(6) 



within the integration volume f2 of each initial conden- 
sate. It is important to emphasize that this argument 
does not rely on the presence of true long range phase 
coherence. In particular, it is valid for 2d Bose gases 
at finite and Id Bose gases at zero temperature, where 
the one particle density matrix approaches a finite value 
no on scales much larger than the interparticle spacing. 
The eventual algebraic decay to zero only appears at dis- 
tances beyond a phase coherence length 1$ that is still 
much larger. For an integration volume that contains a 
large number of particles N 3> 1 and identically prepared 
samples, therefore, the operator a is just the square of 
the condensed atom number in each sample. The eigen- 
values a of a, which are the experimental observables 
according to the standard rules of quantum mechanics, 
thus may take any value between zero and N 2 , where 
N is the number of atoms in either of the two samples. 
The measured integrated density as a function of z varies 
between 2N - 2^/a and 2N + 2yfa. The visibility V is 
therefore simply y/a/N, or 



V 2 



a 

IP 



No 
N 



(7) 



The measured distribution of the visibility thus directly 
reflects that of the condensate fraction. 

Our aim in the following is to calculate the probability 
distribution for the different positive eigenvalues a of a. 
For a given many-body state characterized by a density 
operator p, the associated probability distribution is just 
p(a) — (a\p\a). Mathematically, it is more convenient to 
calculate its characteristic function 



/CO 
dap(a)e laa = Tr [pe la& ] = (e la& ) 
OO 



(8) 



To proceed further, we replace e lcra by its normal-ordered 
counterpart : e %aa : . This approximation neglects commu- 
tator terms which describe the effect of atomic shot noise 
[2"2"] and are of relative order 1 /N [5] . Since typical in- 
tegration volumes for a measurement of the interference 
contrast contain N ~ 10 3 atoms at least, this approxima- 
tion is valid up to corrections of less than a percent. A 
convenient representation for the calculation of the char- 
acteristic function (18]) is obtained by evaluating the trace 



in terms of coherent states. This gives rise to a functional 
integral 

p{a) = (AfjXa)" 1 1 23(^1, Vi)X>(^2,^>) 



x g-SiW-i^il-Safe*] exp 



d d ripi (r)ip 2 (r) 



(9) 



over bosonic c-number fields tp{r^r), which are periodic 
in the interval r = [0, /3] (Here and in the following, we 
adopt units in which h = fes = 1)- Here Si^ are the re- 
spective actions for the interacting Bose gases 1, 2, while 
Mi — j T>(ipi,ipi) exp{—Si[ipi,ipi]) are normalization fac- 
tors. Due to our normal-ordering approximation, the 
fields in the last exponential do not vary with r, in con- 
trast to the fields appearing in the action, but are evalu- 
ated at r = 0. For notational simplicity, ipi(r, 0) = ipi(r). 
By a simple redefinition of a cr/N 2 , Eq. ([9| gives the 
characteristic function for the square V 2 of the visibility 
which is a direct measure of the interference contrast. 

Within the functional integral, it is convenient to 
switch to the density-phase representation ipi(r) = 
\J n i {r)e lipi ^ for the c-number fields. In this represen- 
tation, the square of the visibility V 2 reads 



V 2 = j d d rd d r' Mr)Mr')Mr)Mr') 

^ j d rf rdV v /ni(r)n 1 (r>2(r)ra 2 (r') 

< e i{l<f>2{r)-tpx{r)]-[<p 2 (r')-tpx{r')]} _ Mq\ 



In particular, for temperatures low enough that the influ- 
ence of density fluctuations around an average n{r) may 
be neglected, the visibility 



j d d rn{r)e t[ip2{r) - ipi{r)] 



(11) 



only depends on the phase difference cj> = tf2 — Vi- 



lli. INTERFERENCE CONTRAST AT LOW 
TEMPERATURE 

In the following, we will assume that the two gases 
are identical and describe each one using the quantum 
hydrodynamic action 

n s( r )rv7 . , . r> 



S[ W ]-jfdr/d-r{^[V W (r 



2.9 



[<9 T <^(r,T)]< 



(12) 



Here (3 = 1 /T is the inverse temperature, m is the atomic 
mass and g is a coupling constant, which is just the in- 
verse of the compressibility k. Moreover, n s is the super- 
fluid density, which is inhomogeneous in trapped gases. 
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The action ( 12 ) provides a completely general low-energy 



description of superfiuid Bose gases. In particular, it de- 
scribes 3d gases below the critical temperature for Bose- 
Einstein condensation, 2d gases below the Berezinskii- 
Kosterlitz-Thouless transition [TU [TJ] and also Id gases 
at zero temperature. 

Since V 2 depends only on the phase difference, it is 
advantageous to switch to a new set of variables: 



$ = 



ip 2 + ip x 



ip 2 - ip 1 



(13) 



In terms of these variables, the total action can be rear- 
ranged as 



Substituting the expansion (17) into the action yields 
a diagonal quadratic form 



S 
2 



A^O y v 7 



Introducing the dimensionless variables 



*A = s\\ —tank 



(fax\ 



(18) 



V 2 / 

the characteristic function can then be rewritten as 



(19) 



S = S[cp 1 ]+S[<p 2 ] =2S[$\ + -S[ 



(14) 



which implies that the contribution depending on the to- 
tal average phase $ cancels out in equation ([9]). The 
characteristic function ^ can then be written as 



T>4>e 



-S[<t>}/2 



exp 



N 



n(r)e l 



( 15 ) 

with Af = J T>4> exp(— S[(j>}/2). Note that a constant con- 
tribution in <j> has no effect on the result since it cancels 
out when taking the modulus square. In the following, 
we will thus restrict our analysis to functions without 
constant component, i.e., J d d r <fi(r) = 0. 

For explicit calculations, we follow the technique used 
by Imambekov et al. [HIS] and parametrize the functional 
integral ( 15 1 by expanding <fi in terms of the solutions of 



the imaginary time Euler-Lagrange equation associated 
with the action (12), 



01 



■ [n s (r)W(j)(r,T)] 

m 







(16) 



supplemented with appropriate spatial boundary condi- 
tions. The bosonic nature of the field <p{r,r) requires 
that the solutions satisfy <p(r, 0) = (j)(r, f3). A separation 
ansatz readily gives a family of solutions ij)\{r)e *- T ', 
where A is a formal index labelling the eigenmodes. To 
satisfy the boundary condition on r, we use the expan- 
sion 



A^O 



;/3w A -|- 1 e /3w A _|_ ]_ 



(17) 



Factors have been chosen so that the parentheses evaluate 
to unity at r = and r = (3. 

Unless excluded by the boundary conditions, the 
Euler-Lagrange equation always permits a solution ipo{r) 
which is constant and non-zero in space. Since 
the modes ip\ form a complete orthonormal system, 
J d d r ip\(r)ip\i(r) = S\,\', this implies that all further 
modes have a vanishing spatial average. The condition 
J d d r <f> = therefore translates to the omission of the 
constant mode A = 0. 



= n 

A^O 

where 



dt A e ^ 
x/2^ 



exp ia 



d d r 

IsT 



n{r)e lh ^^ {r) 



(20) 



h {tx} {r) = ct>{r,Q) 




( fax 



coth , 
lux V 2 



t\4>\{r) 



(21) 

is the parametrized iso-phase surface. This expression 
has already been derived by Imambekov et al. 0, start- 
ing from the relationship between the moments of V 2 and 
higher-order correlation functions. The fluctuating sur- 
face h emerged as an abstraction in their paper. Here, 
we see that it is just the shape of the iso-phase surfaces. 
Eqs. ( 20 ) (or rather, the expression for V 2 appearing 



therein) and (21 ) are well suited for numerical use: to ob- 
tain a given realization of V 2 , one generates a large num- 
ber of Gaussian deviates for the amplitudes t\, constructs 
the corresponding surface h and numerically calculates 
the integral to obtain V 2 . Repeating this sequence with 
different sets of amplitudes yields histograms for the pos- 
sible values of V 2 , whose shape approaches the actual 
probability distribution as the number of iterations grows 
large. 

Analytical results can be obtained in the limit of high 
contrasts, 1 - V 2 « 1, where the iso-phase surfaces are 
smooth and one may expand the exponential inside the 
definition of V 2 . As is clear from equation (21), the ex- 
pansion parameter is given by 



max \ / — coth 

A,r V UJ\ 



fax 

2 



(22) 



In all cases that we will encounter in the following, this 
maximum is found for ui m i n = min^o lux ■ 

In the following, we will focus on the regime where 
the reduction of the visibility is dominated by thermal 
fluctuations so that quantum fluctuations (caused by in- 
teractions) may be neglected. Considering, in particular, 
a 2d Bose gas whose characteristic size l z in the direction 
of transverse confinement obeys l z 3> a s , the effect of 
zero point fluctuations of the phase on the reduction of 
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the visibility can be determined by a Bogoliubov calcu- 
lation, which gives [21] 



(V 2 )(T = 0) = 



«o(0) 



= 1 - 



52 
2tt 



(23) 



Here, g 2 = fng<i = \%ira s /l z is the dimensionless 2d 
interaction constant, which has typical values g = 0.1 
[2 117) . which implies practically unit visibility at zero 
temperature. By contrast, at finite temperature, the in- 
tegrated visibility 



(V 2 )(T)-(V 2 >(0) 



l-27j(T)log 



(24) 



decreases logarithmically with the size of the system, 
which reflects the absence of long range order at finite 
temperature in the thermodynamic limit 23, 24J. Since 
the size dependence of the thermal depletion is only log- 
arithmic, finite condensate fractions may be found at low 
enough temperature for realistic system sizes. 

In actual experiments, the temperature is usually 
large compared to the typical frequencies u\ which 
are of the order of the chemical potential fj,. Then 
coth(/3wA/2) can be replaced by 2T/uj\ and hence e 2 = 
2gTmax\ r ip \{r) 2 / ui\. Within this approximation the 
effect of thermal phase fluctuations scales linearly with 
temperature. For weakly interacting 2d Bose gases, 
the necessary condition T ^> /i is well obeyed even in 
the deeply degenerate regime because fi/T — g2/2n at 
= 0(1) [5J. Note that for g 2 < 1 the conditions 
T [i and e 2 <C 1 are simultaneously satisfied for the 
typical phase space densities nA^ that are reached in 2d 
Bose gases [2"5] . 

Quite generally whether or not e is small, i.e., whether 
a physical regime which permits such an expansion exists, 
depends crucially on temperature and the spatial dimen- 
sion d. Because of normalization, the eigenfunctions tp\ 
scale with the characteristic size L = fl 1 ^ of the system 
as L~ d / 2 , while the eigenfrequencies u>\ will scale as L^ 1 , 
independent of the dimensionality. It follows that e scales 
as L v ~ d l 2 . For concreteness, in a homogeneous system 
with periodic boundary conditions, 



1-d 



mT 



(25) 



For the case of 2d Bose gases, which is the main focus of 
our work, e 2 = 2i](T) /it is independent of system size and 
is determined by the exponent tj{T) = (ngAy) -1 which 
gives the decay of the one-body density matrix. 

For e <C 1, the exponent in equation (20) can be ex- 
panded in the form 



exp 




T,^-l E 



e 2 t x t x , (26) 



A/0 



A,A'^0 



since the variables t\ are of order one due to the Gaussian 
weight factors e~* A / 2 . Within this approximation, an ex- 
act calculation of the distribution functions is possible. 



The inclusion of the terms quadratic in e leads to non- 
trivial distributions instead of the Delta functions that 
result in leading order in e [5] . In the general case of an in- 
homogeneous system with a spatially varying superfluid 
density n s (r), thermal phase fluctuations lead to a re- 
duction of the visibility from unity (or — more precisely — 
from its value at zero temperature) of the form 



A^O ^ Wo ^ 



where is the integration volume and 



'A, A' 



Q 

JV 



(27) 



(28) 



Jx = I d d rn(r)tpx(r) 



are dimensionless numbers. In a homogeneous system, 
I\,\> = <5a, A', and J\ = for any A. In inhomoge- 
neous systems, there may be finite "off-diagonal" values 
for /a, A' and finite values for J\. 

For a discussion of some general features of the statis- 
tics of the interference contrast like the dependence on 
dimensionality and the related issue of self-averaging, we 
focus on homogeneous systems (we will discuss the exper- 
imentally relevant trapped 2d system in section IV C ) . It 
is then convenient to define 



^ 2 > = £s 



A^O A 



(29) 



where e has been defined in equation ( 25 1 and uj\ = 
v\/w m i n . Note that the scaling factor 2/e T ~between 1 — 
V 2 and u is large compared to one. While the visibility 
takes values on the interval [0, 1], the auxiliary variable 
u has values on the interval [0, oo] due to our expansion 
of e l ^ r \ The characteristic function 



A^O 



(30) 

of the probability distribution q(u) for the rescaled devia- 
tion u of the visibility from unity is now readily evaluated 
to be 



?m = n 



i 



(31) 



This evaluation ceases to be straightforward when /a, A' 
is not diagonal, since it amounts to the calculation of the 
determinant of an infinite matrix with a non-trivial entry 
structure. 

The logarithm 



, , , 1 {2ia) s ~ (ia) s . s 

s—1 s—1 



(32) 



G 



of the characteristic function q(a), which is the gen- 
erating function of the cumulants of u, can be ex- 
pressed in terms of the spectral zeta function C{A}( S ) = 
Ea/o(^a) _S °f the eigenfrequencies of the quantum hy- 



divergent at s = 1 in two and three dimensions. Since 



drodynamic action (12 1. In particular, it determines all 
cumulants of the random variable u via 



{u s ) c = r-\ s -i)\a {x} {s) 



Note that this calculation does not depend on the explicit 
form of the eigenvalues and eigenfunctions: the geome- 
try of the system is completely contained in the factor 
between u and 1 — V 2 and the spectral zeta function (for 
systems with diagonal Za,A')- 

The precise form of the spectral zeta function obvi- 
ously depends on the geometry of the system and the 
spectrum that follows from it, but the following prop- 
erties are valid for any homogeneous system: C{a}( s ) is 
a monotonically decreasing function of its argument and 
has a lower bound (which is reached in the limit s — > oo) 
equal to the degeneracy of w m i n (note that this essen- 
tial property cannot be reproduced when one replaces 
the sum in C{a}(s) by an integral). It follows that for 
increasing s, the number of frequencies that make a non- 
negligible contribution to the value of C{A}( S ) decreases 
so that higher-order cumulants will essentially depend 
only on a small number of low frequencies. However, we 
will find that it diverges for s = 1 in d > 2 and must be 
rendered finite by the introduction of a UV cutoff. 



Substituting back from equation ( 29 ) , we obtain the 
cumulants of V 2 — 1 (which are identical to the cumulants 
of V 2 except for the expectation): 



((V 2 - l) s ) c = 



(34) 



As we will see, the spectral zeta function remains finite 
for all s > 2 in all relevant cases, thus determining the 
finite size scaling behavior of all higher cumulants: in d 
dimensions, the sth cumulant scales as L s ( 2 ~ d '. 

Specifically, we consider a homogeneous system in d 
dimensions in a hypercubic volume O = L d with pe- 
riodic boundary conditions. Then, LOy~ — c|fe|, with 
k = (2tt/L)(1 1 , . . . where li € Z, and the speed of 
sound c = yj gn s /m. The resulting spectral zeta function 
then reads 



C{k}(s) 



1 



(If 



n=l 



A d {n) 



1 



(35) 



where the prime on the sum indicates that the point 
l\ = ■ ■ ■ = Zrf = is omitted and Ad(n) is the number 
of possibilities to represent the integer n as a sum of d 
squares (including squares of negative numbers). The 
representation on the right-hand side is a special case 
of a Dirichlet series which arises in connections between 
number-theory and modular forms |26) . 



The expectation is then given by equations (34) and 
1. However, simple power counting re- 



( 35 1 with s 



the quantum hydrodynamic action ( 12 1 is an effective 



low-energy description, however, this divergence is an 
artifact. It can be avoided by introducing a cutoff at 
a maximum momentum A = 27r/£, where £ is the healing 
length. The necessity of an explicit UV cutoff has the im- 
portant consequence that the expectation does not follow 



(33) the scaling with system size announced in equation ( 34 1 



in 2d, the cutoff introduces a logarithmic dependence on 
system size that would otherwise be absent so that the 
expectation takes on a non-universal character. During 
the remainder of this article, we will focus on the higher 
cumulants (from the variance on) which are universal. 

Since the variance is independent of system size (as 
are all higher cumulants) while the expectation vanishes 
logarithmically as L — > oo, fluctuations are not self- 
averaging in 2d and the regime 1 — V 2 <C 1 can only be 
reached in finite-size systems (even if they may in fact be 
quite large due to the weak logarithmic size dependence 
of the expectation) . This is in contrast to the 3d situation 
where as a consequence of true long range order, the ex- 
pectation is finite in the thermodynamic limit while all 
higher cumulants decrease with increasing system size. 
In this case, the visibility is a self-averaging observable, 
i.e., for large integration volumes fl, the value obtained 
in a single run is equal to an average over many runs. 
In section |IV| we will focus on system sizes where the 
condensate depletion remains small and the visibility is 
close to unity, i.e., L -C 1^ — ^e 1 ' 217 ^ and the system 
is a true condensate rather than a quasi-condensate [21] . 
The opposite case will be discussed in section |V} 

The fact that all cumulants from the variance on are 



finite and obey the scaling given in equation ( 34 1 implies 
that the probability distribution p(V 2 ) is universal and 
non-Gaussian. The following section is devoted to the ex- 
plicit analytic calculation of this distribution in different 
geometries. 



IV. ANALYTICAL RESULTS IN 2D 

For a given geometry and boundary conditions, one 
can always evaluate the spectral zeta function numeri- 
cally. Equation (34 1 then permits the calculation of an 



veals that the spectral zeta function is ultraviolet (UV) 



arbitrary number of cumulants. However, this is not suffi- 
cient to obtain an analytical expression for the underlying 
probability distribution p(V 2 ) which requires a closed- 
form expression for C{A}( S )- This may be obtained in the 
2d case for simple geometries, which we will discuss in 
this section. 

In order to clarify the applicability of our results to ex- 
periment, it is necessary to define what we mean by "2d" 
in practice. We consider the two gases to be strongly 
confined along the z direction with a trapping potential 
uj z which is sufficiently strong to ensure uj z 3> T and 
w 2 > ;i so that the gases reside in the harmonic oscil- 
lator ground state along the z direction. However, the 
spatial extension l z = y/l/muj z shall still be large com- 
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pared to the scattering length a s . This regime, which is 
sometimes referred to as "quasi-2d" [27], accurately de- 
scribes the situation in typical experiments on cold atoms 
where the tight confinement is realized using an optical 
dipole potential [2J [TSJ [2"8"H3"T] . It is particularly simple 
in that the interaction may be described by a simple di- 
mensionless constant g = mg = \/8wa s /l z (for not too 
strong interactions. Here and in the following, we drop 
the subscript "2" for ease of notation) whereas the inter- 
action constant in a two-dimensional system with l x < a s 
depends on the chemical potential and thus effectively on 
the spatial density [27] . 

The first case we discuss is that of a rectangle with 
periodic boundary conditions. The latter are certainly 
artificial, but the results are important from a concep- 
tual point of view because the calculation can be carried 
out in closed form. We will discuss two limiting cases of 
this particular case (isotropic and strongly anisotropic) 
before moving on to the experimentally relevant case of 
a harmonically trapped sample. 



P(s) 1 0.915966 1 0.968946 1 0.988945 1 0.996158 1 

TABLE I: The first four relevant values of the Dirichlet beta 
function 



where the prime on the sum signifies that the value 
1 1 = I2 = must be omitted (for a derivation of equa- 
tion (38), see [35]). Here, ((s) = J2kLi k~ s ^ s the 
mann zeta function and /3(s) = X^o( — 1)^/(2Z + l) s is 



the Dirichlet beta function. As already stated, (38) is 
undefined for s = 1 since the Riemann zeta function di- 
verges, i.e., the expectation of u cannot be calculated 
without a UV cutoff. For s > 2, i.e., for all cumulants 



except the expectation, ( 38 ) is finite. The introduction 



of a high-fc cutoff will only cause small deviations since 
the sums are infrared-dominated and we may write 



logg(er) = ia(u) +2^ 



(2iaY 



C(s)/3(s) 



(39) 



s=2 



A. Homogeneous square sample 

The conceptually and mathematically most elementary 
case is the one where the gas is homogeneous and confined 
to a rectangular area of extension L x aL, < a < 1, 
with periodic boundary conditions. First, we will focus 
on the particular case a = 1, i.e., a square sample, but 
it is convenient to introduce the notation directly in its 
slightly more general form for arbitrary a. In this case, 
the eigenfunctions take the simple form 



cos(fc-r) 



sin(fc-r) (36) 



,-(l-V 2 ). (37) 



aL 2 y ' ' ^ k V aL 2 

with k = (2tt/L) x (nx/ a,n%), and ni^ € Z. 

We may now give an explicit expression for the expan- 
sion parameter e defined in equation ( 22 ) and the rela- 
tionship between V 2 and the auxiliary variable u. For 
this particular set of eigenfunctions and eigenvalues, 

2 1 rnT 2 2 

e = —2 = —V(T) 

air^ n s ira 

Quite remarkably, e and hence the entire probability dis- 
tribution p(V 2 ) have no explicit dependence on the in- 
teraction constant g, which is hidden in the non-trivial 
relation between the superfluid density that enters 77 (T) 
and the bare 2d density n. For weakly interacting Bose 
gases, this relation has been worked out in [3"2"] . 

As we already stated, the calculation of a closed-form 
expression for the probability distribution p(V 2 ) requires 
a closed- form expression for q£\ ■ If turns out that in two 

dimensions and in the absence of a high-/c cutoff is 
a special case of a lattice sum first calculated by Lorenz 
[33] and Hardy [34] who showed that 

, v — 1 



{fc}> 



(it + qy 



= 4C(s)/3(s) 



(38) 



where (u) is evaluated using a finite cutoff. The different 
UV-behaviour of the expectation and the higher-order 
cumulants has the consequence that the probability dis- 
tribution for u is universal except for a cutoff-dependent 
shift. Thus, the variable u — (u) has values on the entire 
real axis. As long as the resulting probability distribution 
takes on negligible values for values of u corresponding to 
visibilities outside the interval [0,1], this does not create 
any inconsistency. 

As pointed out by Bramwell in the context of the or- 
der parameter distribution for the 2d XY model in the 
low temperature limit [35] , ( 39 1 is closely related to the 
Gumbel distribution 



Pg(x) = exp[-(z + 7) - e ( x+ ^} 



(40) 



(with the Euler constant 7 = 0.577 . . . ) which determines 
the statistics of the interference contrast for weakly inter- 
acting Id Bose gases at zero temperature [7J[5] (note that 
( 40 1 is a normalized distribution with zero average and 



variance tt /6). Its cumulant generating function reads 



z — ' s 

s=2 



(41) 



In fact, there are two non-trivial differences: (i) the pres- 
ence of the Dirichlet beta function in (39). This func- 



tion rapidly converges to one so that it may be replaced 
by unity to a good approximation (see table [I]). Since 
the expectation is non-universal anyway, this essentially 
amounts to increasing the variance by 9%, or equiva- 
lently, the width of the distribution by 4% (this may 
be seen as a convolution with a normalized Gaussian of 
appropriate width), (ii) the global factor of two which 
implies (accepting f3(s) ~ 1) that q(u) is the convolu- 
tion of two identical Gumbel distributions, and p(V 2 ) its 
scaled mirror image. Thus, there is a striking similarity 



between the 2d case at finite temperature discussed here 
and the Id case at vanishing temperature. However, ow- 
ing to the different corresponding eigenspectra, passing 
from the latter to the former case does not amount to 
the simple replacement K i-> 1/2-q as suggested in [5]. 



T 




u — (u) 




(«-<«»M^ 

FIG. 3: The same data as in the previous figure (except for 
the choice of shown temperatures), but with normalised vari- 
ance. In this representation, the convergence towards the 
asymptotic shape is considerably faster. 



FIG. 2: Numerically calculated distributions q(u — («)) for 
T c /T = 1,2,4,8,32 (symbols, cf. legend) and a convolution 
of two Gumbel distributions (continuous line). The lines are 
guides for the eye. Here and in the following plots, the sta- 
tistical error is of the order of the symbol size. 



The evolution of the distribution q(u) with decreasing 
temperature, obtained numerically, is shown in Fig. [2] 
For simplicity, we have subtracted the non-universal ex- 
pectation so that all distributions are centered around 
zero. The numerical data was obtained by generat- 
ing 100 000 random surfaces per curve using a total of 
1000 modes on a 16 000 point grid and then calculat- 
ing V 2 using equations (20) and (21) without any ap- 
proximation beyond the replacement of coth(j3u}\/2) by 



2T/u>\. In particular, there is no expansion of exp(ih) 
to second order in h. For sufficiently low temperatures, 
the distribution approaches the universal low tempera- 
ture distribution, whose characteristic function is given 
in equation ( 39 1 . The numerically calculated distribu- 
tion for n s /mT = 2/tt should of course not be taken seri- 
ously: while the superfluid density remains finite at the 
critical point, the range of momenta where the quantum 
hydrodynamic action (12) is valid approaches zero. A 



proper result for the distribution of the interference con- 
trast near T c requires calculating the full counting statis- 
tics of the condensate number near the BKT-transition, 
as will be discussed in section fVl 

As one can see in Fig.[3j the actual shape of the distri- 
bution that is obtained from scaling q(u) by its standard 
deviation yj (u 2 ) c is already quite close to the asymp- 
totic result for n s /mT = 2/tt and converges rapidly with 
decreasing temperature. 



B. Strongly anisotropic rectangle 

While there is no general closed form expression for 
C{k} f° r arbitrary values of a, it is nonetheless possible to 
obtain an analytical expression for p(V 2 ) in the limiting 
case a -C 1, i.e., for a very anisotropic sample. The rea- 
son lies in the mathematical structure of the cumulants: 
except for the non-universal expectation, all cumulants 
are given by sums which are dominated by the lowest- 
lying eigenvalues. Already for moderate values of the as- 
pect ratio a (around 1/5. Preliminary experimental data 
has been taken for aspect ratios even lower, < 0.1 |17j). 
the contribution of the modes along the shorter direction 
of the samples becomes negligible and one obtains 



2C(2 S ) 



(42) 



It is important to keep in mind that in order to re- 
main in the two-dimensional regime, the aspect ratio 
must not become too small, i.e., equation (42) holds 
under the condition that a -C 1, but at the same time 
h 2 /m(aL) 2 < n,T. If the second inequality is violated, 
the system becomes effectively one-dimensional. At the 
same time, e as given in equation (22) ceases to be a 



small parameter so that one can no longer justify the ap- 
proximate treatment of e l ^ r \ However, we emphasize 
that equation (42) is well satisfied (for s > 2) already for 
moderately small a so that the strongly anisotropic 2d 
regime is well-defined. 
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Substituting ( 42 1 into equation ( 32 ) yields 



ig (it) — 



or (defining ii n 




73) 



(44) 



This function is meromorphic in the entire complex 
plane since the branch cuts of the square roots before 
and inside the cosecant cancel each other. This per- 
mits to explicitly calculate its inverse Fourier transform 
q(u) = (27r) _1 da q(a)e~" La using the residue the- 
orem. In the upper half plane, q(a) has no poles and 
falls off as exp[— ^/2Im(cr)] for large arguments. Thus, 
for u — M m ; n < 0, one may close the integration contour 
with a half-circle over the upper half plane and the in- 
tegral vanishes. For u — u m in > 0, one must close the 
contour in the lower half plane where q(a) has an infinite 
number of poles a n = -in 2 /2. There, e~ la ^ u ~ Umln ^ q(a) 
has the residues R n = (-l)™" 1 ^ 2 e~" 2 / 2 . Thus 
we obtain 



q{u) 





X-i(-i) 



u < u u 

Tl 2 



This can be written in a more compact form in terms of 
its cumulative distribution function: 

q(u) = 6(u (e(— ")/ 2 ) , (46) 

where 9 is the Heaviside function, and $i(z) = 1 + 
2^^ 1 (— l)™z™ is a Jacobi theta function. As a short- 
hand, we will refer to the distribution described by equa- 
tion (45) as "Jacobi distribution". 



The evolution of the probability distribution with de- 
creasing temperature in the strongly anisotropic case is 
shown in Fig. [1] (for a = 1/10). Unlike the isotropic case, 
the value V 2 = is actually quite probable at tempera- 
tures close to the critical point so that the shape of the 
distribution shows a clear qualitative change as the tem- 
perature is lowered towards the asymptotic regime. This 
observation is in qualitative agreement with preliminary 
experimentally data taken at ENS [T7j . Note that as long 
as the probability for a vanishing visibility stays finite, 
the shape of the distribution depends on the expecta- 
tion and is thus explicitly cutoff-dependent. However, 
since this dependence is only logarithmic we expect the 
qualitative evolution of the shape to be insensitive to the 
precise value of the cutoff. 

In turn, Fig. [5] shows the evolution of the probability 
distribution in the asymptotic low temperature regime 
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FIG. 4: Probability distributions p(V ) of the visibility as a 
function of V 2 /{V 2 } for T c /T = 1,2,4,8, for an anisotropy 
a — 1/10. Unlike the distribution in the isotropic case, 
the distribution for a strongly anisotropic sample undergoes 
strong shape modifications as T c /T increases. 




FIG. 5: Probability density q(u) for a — 1,1/2,1/5 as well 
as the asymptotic distributions for the isotropic (heavy line) 
and strongly anisotropic (thin line) regime. For a < 1/10 
(not shown), the results become indistinguishable from the 
strongly anisotropic limiting case. A value of T c /T — 256 has 
been chosen to ensure that the distributions are well within 
the asymptotic low temperature regime. 



(T = T c /256) with changing aspect ratio, confirming our 
earlier statement that the strongly anisotropic regime is 
reached already for a < 1/5. 



C. Harmonically trapped sample 

As we have pointed out in the preceding discussion, the 
cumulants (and hence the shape of the distribution) are 
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dominated by the excitations with the lowest frequency. 
The periodic boundary conditions we have used up to this 
point are thus quite artificial: even for a homogeneous 
system, going over to different boundary conditions will 
have effects on the shape of the distribution (see [5] for 
examples). However, the differences are mere numerical 
factors appearing in the cumulants (for an example, the 
variance in a 3d homogeneous sample is smaller by a 
factor of 1.67 when Dirichlet boundary conditions are 
used instead of periodic ones [37]), the dependence on 
physical parameters remains unaltered. 

We now consider the geometry most relevant for actual 
experiments: a sample which is harmonically trapped. 
For simplicity, we consider an isotropic trap with trap- 
ping frequency uj. If the particle number is large enough 
to warrant Ng ^> 1 [B] (which is readily fulfilled for 
g ~ 0.1 and TV ~ 10 3 ) and the temperature is sufficiently 
low, the density distribution takes the form of a Thomas- 
Fermi profile 



n(r) = n(0) 



1 



TV 



9{R-r) 



R = 



In the following, we will assume the entire sample is su- 
perfluid and we need not distinguish between the super- 
fluid and the total density. This requires the tempera- 
ture to be substantially below Tbkt = (4-kNlj 2 /gD 2 ) 1 / 2 , 
where the phase space density D — nX 2 ^ reaches the 
critical value D c = log(380/g) of the BKT transition 
in the trap center [38 . Moreover, we want to be in 
the regime of near-unity visibility so that we must fulfill 
5(j>{R) = (n(0)A|)- 1 log(i? 2 /C 2 ) < 1 IS, or equivalent^, 



1 



log(2giV/7r) 




(48) 



For N = 5000, ~g = 0.1 and w = 2tt x 20 Hz, one ob- 
tains Tbkt = 95 nK and T$ = 132 nK (for lower in- 
teraction constants, both are even higher) so that the 
low-temperature regime we are considering is within ex- 
perimental reach. 

In analogy to the calculation of the low-energy modes 
in 3d [39 , the solutions Vxz( r j &) 01 the Euler-Lagrange 
equation ( 16 ) with open boundary conditions yield the 
frequencies 



n(n + 2) - P 



(49) 



Here, n = 0, 1, 2 . . . is a radial and / = — n, — n+2, . . . , n— 
2,n an azimuthal index |40j . The eigenmodes are of the 
general form 



ipn,l(r,9) = 



P n ,\i\(r/R) 



R 



COS(W)/y/TT l>0 

1/V2n 1 = , (50) 

sm(\l\e)/y/w l<0 



Efc cl^^x* are polynomials the co- 
efficients of which may be obtained from the recursion 



where P n ,\i\(x) 



relation 

,>,!*!) r;2 



[l 2 - (Jfe + 2) 2 } = a^ W) [n(n + 2) - k(k + 2)] . (51) 



l k+2 



Hence, the P n n\(x) are either even or odd, and the high- 
est and lowest occurring power of x are n and respec- 
tively. The magnitude of the lowest coefficient is fixed by 
the normalization condition J Q dx xP n ^\i\ (x) 2 = 1. 
The eigenmodes satisfy the orthogonality relations 

d6 drrVvM)Vv,Kr,0)=<W<*M' (52) 
o Jo 

on a disk of radius R. For different I, the orthogonality is 
assured by the azimuthal part, for equal I, by the radial 
part of the eigenfunctions. 

A particularly important subset of the eigenfunctions 
is formed by the surface modes 



ipn,±n(r) = \/2(n + l) 



\ cos(n0)/ v /7r I = n 
2gfi<fi) — — v v ■ 'ijn+i| sillM ^ l = -n 

rnco 2 (53) 
(47) with eigenfrequencies uj n> ± n = ^/nu. By inspecting equa- 



tion ( 49 ) , one finds that most of the lowest-lying modes 



are such surface modes. This is physically intuitive, since 



the phase stiffness as given by equation ( 47 ) is lower close 



to the rim so that low-energy excitations should live on 
the boundary of the sample. It is thus to be expected that 
the condensate fraction and hence the interference con- 
trast is dominated by the behavior of the surface modes. 

If we substitute the eigenmodes and eigenfrequencies 
of the surface modes into the definitions of e and u, we 
obtain (noting that ip n ,i(r) 2 /oj 
\l\ = 1 and r = R) 



2 , takes its maximum for 



4 mT 



it n s (0) 



^(1 



V 2 



(54) 



Once again, as in the homogeneous case, there is no ex- 
plicit dependence on the interaction constant. However, 
in the present case, there is an implicit dependence on 
the interaction constant since the latter defines the geom- 
etry of the sample and the density in the center is given 
by n(0) = mNuj 2 /irg. 



Carrying out the expansion (271, we find that in the 



case of harmonic trapping the integrals ( 28 ) become non- 
trivial. One readily finds that I n ,l,n',V — &l,l'In,n' (|^|) 
and J nt i = 5ifiJ n which follows immediately from the az- 
imuthal part of the eigenfunctions. This shows already 
that there is no / that "couples" different surface modes 
(they all differ in I) and no J that involves any surface 
modes. The diagonal elements for the surface modes are 
readily found to be I n ,±n{n) = 2/(n + 2). On closer 
inspection, one finds that the matrices I n ^ n i(\l\) are tridi- 
agonal in the sense that they are non-vanishing only for 
n' = n, n ± 2 (we recall that n and I are either both even 
or both odd), likewise only J 2 = l/v3 is finite. 

Taken together with the physical intuition that the 
statistics should be dominated by the surface modes, 
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Approximation 


<« >e 




exact 


2.813 


1.677 


diagonal 


2.496 


1.580 


surface 


2.160 


1.470 



TABLE II: Variance and standard deviation of u in the har- 
monically trapped case, exact and in two different approxi- 
mations. 



these results suggest that it should be a reasonable ap- 
proximation to disregard J2 and the non-diagonal ele- 
ments of the matrices I n ,n'{\l\) which would permit to 
carry out the integration leading to equation (31 1. In 



order to see whether this leads to correct results we eval- 
uate the variance of u first exactly, then in the "diagonal" 
approximation. A straightforward calculation gives 



E 



'In,n>{\l\f 



6 z^ 



,(0) 



1 

72 



(55) 



which is reminiscent of a structurally similar expression 
for the variance of condensate fluctuations in a 3d har- 
monic trap by Giorgini et al. |41| . but with two addi- 
tional terms which come from the fmiteness of J 2 . Once 
again, the primes on the sums are reminders that the 
sums go only over permitted values of the indices. 

In table [Til we give the numerical value of equation 



(55) calculated in three different manners: first exactly, 
i.e., without approximations apart from the numeric cal- 
culation, and in two different approximations. In the "di- 
agonal" approximation, we disregard the two last terms 
which are generated by J 2 and take into account only 
the diagonal elements of I n .n'(\l\) in the first term. The 
"surface" approximation goes even one step further by 
dropping all contributions except those coming from the 
surface modes. We see that the diagonal approximation 
is quite satisfactory (we recall that the higher cumulants 
are increasingly dominated by the lowest-lying modes so 
that we expect the approximation to improve with in- 
creasing cumulant order) and even the elementary sur- 
face approximation fares reasonably well. 
In the diagonal approximation, u reads 



r, 



,1 



I W n,l 



(56) 



The approximation on the right-hand side is better than 
within one percent for all s > 2. The global factor of 2 
comes from the two-fold degeneracy of the surface modes. 
As in the strongly anisotropic rectangle case, it compen- 
sates the global factor of 1/2 in equation (32). 

While it does not seem possible to derive a closed-form 
expression for the associated probability distribution, the 
form of the spectral zeta function suggests that the dis- 
tribution should be something intermediate between a 
Gumbel distribution [where the spectral zeta function is 
proportional to C( s )] an( i a Jacobi distribution [where 
the spectral zeta function is proportional to £(2s)], i.e., 
more asymmetric than the former, but less asymmetric 
than the latter. Of course, the contributions from the 
neglected modes will render the actual distribution some- 
what more symmetric than this argument suggests, but 
it should remain qualitatively valid. This is supported by 
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FIG. 6: Scaled distribution q(u) for the harmonically trapped 
case, for T c /T = 1,4,64. A Gumbel distribution (not a con- 
volution of Gumbel functions) scaled to have unit variance is 
shown for comparison (continuous line). The sharp fall-off at 
the right end of the T — T c curve corresponds to zero visibility 
while the T = T c /64 curve is already well within the asymp- 
totic low-temperature regime. Quite amusingly, the profile 
for T — T c /T = 16 (not shown) comes out almost superposed 
with a Gumbel distribution while the asymptotic distribution 
comes quite close to it, but remains more sharply peaked, in 
agreement with the arguments presented in the text. 



and the spectral zeta function is 



n.l V 



(57) 



In the surface approximation, this can be written explic- 
itly as 



-J \n(n + 2) 



H^) C(3s/2) 



(58) 



numerical results which are calculated without approxi- 
mations, as can be seen in figure [6j 



V. DISTRIBUTION AT THE CRITICAL POINT 

A quite interesting aspect associated with the statistics 
of interference amplitudes is the possibility to measure 
the universal probability distribution of the order param- 
eter near a critical point. Indeed, as has been shown in 
section ITT) the visibility is identical with the condensate 
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fraction provided the integration length is much larger 
than the interparticle spacing. Our calculation of the re- 
sulting visibility distribution in the previous section is 
valid deep on the Bose-condensed regime, where the vis- 
ibility is close to one. 

In the following, we want to discuss the situation close 
to the critical point of Bose-Einstein-condensation, that 
is in a regime where the system size L is large com- 
pared to microscopic lengths but of the same order or 
smaller than the correlation length £ of the infinite sys- 
tem. Mathematically, this may be expressed by a dimen- 
sionless parameter 



= tL 1 '" = ± 



l/u 



= 0(1) 



(59) 



that measures the deviation from the critical point due 
to the finite system size. Here t — (T — T c ) /T c is the 
dimensionless distance from the bulk critical temperature 
T c and v « 0.672 the critical exponent that characterizes 
the divergence £ ~ \t\ ~ v of the correlation length in a 3d 
BEC [321 133] • This exponent has in fact been measured 
also in dilute ultracold gases [33] ■ Quite generally, finite 
size scaling predicts that the probability distribution of 
a two-component order parameter s (for a BEC, s 2 = uq 
is the condensate fraction) in the critical regime has a 
scaling form [131 131] 



Here. 



y(d) = (d-2 + V )/2 



(60) 



(61) 



is related to the standard anomalous dimension rj of the 
XY-model while p%{z, x) is a universal, non-Gaussian dis- 
tribution that only depends on z 2 = n L 2y . The exis- 
tence of such a distribution for properly scaled block-spin 
variables sL v with finite moments of arbitrary order, is 
in fact a basic assumption of the renormalization group 
approach to critical points, as emphasized by Parisi [15] . 
Precisely at the critical point, where x = 0, the distri- 
bution is determined by the effective potential V* s (z) of 
the underlying field theory at the fix-point via 



p2(z, a: = 0)~exp [-V* ff (z)] 



(62) 



Within the standard two-component $ 4 -theory, the dis- 
tribution pz(z,x) in the 3d case has been calculated by 
Chen et al. [47], taking into account the singularities 
associated with the Goldstone mode. The resulting dis- 
tribution at the critical point x = is shown in Fig. [7] 
It exhibits a maximum at z 2 ~ 1.146. The most prob- 
able value for the number of particles in the condensate 
therefore scales as N ~ L 2 ~ v which is sub-extensive, as 
expected right at T c . The average typical visibility 



(63) 



at the critical point will therefore vanish with an anoma- 
lous power of the integration length L, i.e., basically like 



1/L because r\ ~ 0.03 is rather small. Moreover, the fact 
that the variance of the variable z 2 is a universal con- 
stant of order one, the so-called Binder cumulant |45j . 
implies that the fluctuations of the visibility scale with 
the same anomalous power of the integration length as 
the average. 

In the 2d case, there is no breaking of a continu- 
ous symmetry at finite temperature. Instead, there is a 
Berezinskii-Kosterlitz-Thouless transition to quasi long 
range order below T c , where the gas is a proper super- 
fluid but not a BEC. The absence of a finite correlation 
length below T c in this case does not allow to define a 
simple analog of the variable x in ( 59 1 . Right at the 



BKT-transition, however, one expects again a universal 
order parameter distribution function p^C 2 ) f° r the vari- 
able z 2 = no-L 17 , since y = rj/2 in two dimensions. In 
contrast to the situation discussed in section |IV[ where 
the distribution of the interference contrast has been cal- 
culated deep in the superfluid regime and the visibility 
is close to one, the distribution p%( z ) with its anomalous 
scaling applies to 2d Bose gases whose size is much larger 
than the phase coherence length l^,. The thermal phase 
fluctuations then imply an average condensate fraction 
(no) ~ i -1 * which decreases with system size. The typi- 
cal value yj (V 2 ) ~ L^ 11 of the visibility is therefore close 
to zero. In fact, since the 2d superfluid phase corresponds 
to a line of critical points at any T < T c , this behavior 
of the average visibility is valid at arbitrary temperature 
below T c in the limit L — ¥ oo with a temperature de- 
pendent exponent Tj(T) which reaches its critical value 
rjc = 1/4 at T c . Note that, independent of the precise 
form of the distribution p^iz), the very existence of a 
scaling variable z 2 = n^L 11 immediately implies the sub- 
extensive scaling (No) ~ L 2 ~' n of the average number 
of particles in the condensate for an interacting 2d Bosc 
gas and its anomalous fluctuations VariVo ~ (Nq) 2 [35] 
in the thermodynamic limit. 

In order to observe this anomalous scaling, the sys- 
tem size must be large compared to the phase coherence 
length, L > I $ = ^e 1/2ri(T \ At the critical point, this 
is readily fulfilled for typical system sizes of the order of 
some 10um and healing lengths of the order of 0.1um. 
By contrast, for T <C T c , the system size required to be 
in the anomalous scaling regime rapidly exceeds exper- 
imentally feasible values. Therefore one has L <C i<j> in 
practice and the visibility distribution can be deter min ed 
by an expansion around V 2 ~ 1 as done in section IV 



VI. CONCLUSION 

In conclusion, we have shown that interference ex- 
periments may be used as a direct measurement of the 
statistics of the condensate fraction in ultracold Bose 
gases. Unlike in interference experiments in classical op- 
tics, where the fringe visibility is determined by a de- 
terministic cross correlation function of the optical fields 
[4Tfl [5U] . the interference contrast of matter waves is a 
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FIG. 7: The universal distribution ^3(2,0) for the order pa- 
rameter z in 3d directly at the critical point as calculated in 

12. 



quantum observable. Repeated experiments with identi- 
cally prepared condensates therefore produce a statisti- 
cal distribution of values instead of a reproducible single 
value. The resulting distributions are non-Gaussian even 
in the thermodynamic limit and have been calculated ex- 
plicitly for 2d Bose gases at temperatures such that their 
effective condensate fraction is close to one. Quite gener- 
ally, the interference contrast is a self-averaging observ- 
able in situations with long range phase coherence. Our 
findings for the 2d strongly anisotropic case are in qual- 
itative agreement with preliminary data taken at ENS 
[17) . Clearly, a quantitative comparison between theory 
and experiment is needed to verify our predictions. In 
particular, the interference statistics might be used as 
a precise thermometer of the gases, similar to what has 
been achieved in Id gases [TO]. A quite interesting open 
problem, both from a theoretical and an experimental 
point of view, is the analysis of the interference contrast 
near the transition to the normal phase. It offers the pos- 
sibility to directly measure the distribution of the order 
parameter 1 near the critical point, a quantity that is very 
hard to measure otherwise. 
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Appendix A: Universal scaling in 3d 

As was shown in section [TTJ the distribution of the con- 
densate number, which is related to the intensive two- 
component vector order parameter s that describes Bosc- 
Einstein-condensation from the point of view of statisti- 
cal physics by N = L 3 s 2 is not a simple Gaussian. This 
is a result of the fact that in the case where the broken 
symmetry is continuous, the order parameter correlation 
length is infinite for all temperatures below T c |37j . The 
universal distribution function of the condensate number 
below T c is in fact contained in the result (32) for the 



logarithm of the characteristic function of the random 
variable u = 2(1 — V 2 )/e 2 . In 3d, the small parameter e 



defined in equation (25 1 can be written in the form 



1 frcn 

TT 2 L 



(Al) 



where we have introduced the Josephson length £j = 
mT/n s . For any finite temperature therefore, e goes to 
zero for a system size L much larger than the Josephson 
length. The universal distribution p{u) for the fluctuat- 
ing variable 



2tt 2 L 



1 



Nl 
N 2 



(A2) 



which determines the distribution of the condensate num- 
ber of a 3d BEC below T c is fixed by the exact cumulants 
given in ( 33 ) . It depends on the 3d spectral zeta function 



C{A}(*)= E 



1 



(if + q + qy 



J2 Mn) 

n=l 



II " 



(A3) 



for which, unfortunately, no closed-form expression seems 
to exist [35]. It is evident, however, that C{a}(s) is con- 
vergent for all s > 3/2 and thus all cumulants except the 
first are finite. The variable with a proper, non-Gaussian 
distribution in the limit L — > 00 is thus n 2 , • L/£j which 
implies that the condensate fraction is a self-averaging 
variable. Its fluctuations, however, are not of order l/Q 
as usual but only decay like 1/L 2 [37]. Since all higher 
cumulants including the variance are finite and have the 
same scaling with system size L, the ratios (u s ) c / {u 2 ) s J 2 
are constant and finite. A special case of this result has 
in fact been found by Kocharovsky et al. 51 , who calcu- 
lated the cumulants of the number of condensed atoms in 
a 3d BEC within a Bogoliubov approach. In particular, 
to leading order in e, our cumulants from equation (34) 



agree with theirs, showing the close connection between 
the condensed fraction and the interference amplitude. 



1 Note that for 2d gases there is no true order parameter, yet 
there is a nontrivial distribution of the number of particles at 
zero momentum. 



Appendix B: The Id case at zero temperature 

In this appendix, we discuss the case of a homoge- 
neous Id Bose gas at vanishing temperature. Using 
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c = y/ ' gn s /m and substituting n s /m = cK/tt, where 
K is the dimensionless Luttinger parameter, the hy- 
drodynamic action ( 12 ) governing the phase difference 
(\> = ip 2 — fi of two interfering Id Bose gases has the 
form 



So[0] = 



K 
Aire 



dx J dr [[d T ^] 2 + c 2 [d x <j,} 2 j , (Bl) 



where we have kept f3 finite. Using the Fourier expansion 

1 oo 

<j>{ X ,r) = E 4>kb>n)j {hX - UnT) > (B2) 



k ri= — oo 

where k = 2irl/L with Z g Z, and w„ = 2im/f5 are the 
bosonic Matsubara frequencies, the action takes the di- 
agonal form 



K 

4-7TC 



E(c 2 fc 2 +^)|^K)| 2 . (B3) 



The generating function p(a) 



for the square 



of the visibility requires calculating a functional integral 
with a perturbation 



Si = j^2 I dx I dx ' n{x)n{x') cos [<j>(x) - <f>(x')] (B4) 



to the action (Bl I. This perturbation only contains the 
phase difference 4>{x) = cj>(x, 0) on the boundary in imag- 
inary time t. Except for 4>{x) = 4>{x, r = 0) all variables 
are therefore Gaussian and can be integrated out. The 
problem then is completely analogous to that of backseat - 
tering from a single impurity in a Luttinger liquid dis- 
cussed by Kane and Fisher [52"] . 

Upon elimination of the modes 4>k (r ^ 0), one obtains 
the reduced free action 



SM = 



K 
2k 



(B5) 



fc=-A 



for the remaining, non-Gaussian degrees of freedom, 
where we have explicitly written the ultraviolet cutoff 
A. This corresponds to a non-local action in space of the 
form 



Sq[<Kx)] = 



8tt 2 



dx / da;' 



(x) ~ 0(1') 



(B6) 



that arises in r space for dissipative quantum mechan- 
ics of a single particle |53j or in the study of nontriv- 
ial ground states of open strings [53] . Equations (B5) 



and (B6| are equivalent if the associated dimensionless 



strength a of the dissipation is related to the Luttinger 
parameter by a = 2K . 



Following the arguments of section IV the distribution 
of the interference contrast can be calculated analytically 
in the limit e 2 = 1/K <C 1 by expanding S% in equation 



(|B4|) to second order in 0. Again, it is then natural to 

d «t(1-V 2 )\ 



consider the characteristic function q{u) = {e" 
which corresponds to a perturbation (for a homogeneous 
system with n(x) = N/L) 



SM = - 



2L 2 



dx / dx' [0(x) - <j)(x')Y 



(B7) 



Substituting the Fourier series representation for <p(x), 
this becomes 



L 



(B8) 



Mo 



The functional integral is now Gaussian and can be eval- 
uated exactly, giving 



*>=n(i-fsr- 



fe>0 



(B9) 



logg(a) = -(l-V 2 



E 

s=2 



1 / ia 
K 



(BIO) 



where again the expectation is explicitly cutoff- 
dependent. Comparison with equation (40) shows that 
the variable K{(V 2 ) — V 2 ) has a Gumbel distribution of 



the normalized form given in equation ( 40 ) as derived in 



Note that the action Sq + Si as given by equations 



(Bl) and (B4| differs from the action of the boundary 



sine-Gordon model that appears for dissipative quantum 
mechanics in a (purely imaginary) periodic potential. In- 
stead it corresponds to a classical Id XY-model model 
with infinite range interactions. However, a mapping to a 
sine-Gordon model (relying on a Hankel transform rather 
than a Fourier transform of the probability distribution) 
is possible in the thermodynamic limit and has been used 
by Gritsev et al. in [7] to calculate the distribution func- 
tion of the interference contrast for arbitrary values of 
the Luttinger parameter K. 
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